## replication script for: ##
## 'Allocating Costly Influence in Legislatures' ##
## Provins, Monroe, and Fortunato ##
## Journal of Politics ##
## R version: 4.0.2 ##
## date of script construction: June, 21, 2021 ##

## load libraries ##
	library(stargazer)
	library(mnormt)
	library(lme4)
	set.seed(0624)

## load data (remember to set your directory!) ##

 	setwd('')
 	data <- read.table('finalData.txt', sep = '|')
	summary(data)
  
## minorityCommitteeShare = share of committee seats allocated to minority ##
## minorityChamberShare = miority's total chamer seat share
## proportionalityRule = indicator for proportionality rule
## minorityLeaderConference = indicator for rule providing minority leader a formal consultation
## majorityCohesiveness = standard deviation of majority ideal points * -1
## proceduralPower = factor analytic summation of majority's procedural power
## totalCommitteeSeats = logged number of total seats on committee
## unitId = identifies legislature, also identifies session for congressional observations
## committeeName = the name of the committee

## theory picture (Figure 1) ##

	x <- seq(0, 0.5, by = 0.01)
		
pdf('newTheoryPic.pdf', width = 10, height = 10)
	par(mfrow = c(2, 2))
	plot(1, 1, type = 'n',
		xlim = c(0, 0.5),
		ylim = c(0, 0.5),
		xlab = 'Minority Chamber Seat Share',
		ylab = 'Minority Committee Share',
		main = '',
		axes = FALSE)
	axis(1, at = seq(0, 0.5, 0.1))
	axis(2)
	y <- x
	lines(x, y, lwd = 2, lty = 2)
	lines(x, y, lwd = 6, col = rgb(0, 0, 0, 0.5))
	legend('bottomright', legend = c('Predicted Minority Share'), lwd = 6, col = rgb(0, 0, 0, 0.5), bty = 'n')
	text(0.05, 0.45, substitute(paste(alpha, ' = 0')))
	text(0.05, 0.4, substitute(paste(beta, ' = 1')))
	mtext('Proportional Allocation', font = 2)

	plot(1, 1, type = 'n',
		xlim = c(0, 0.5),
		ylim = c(0, 0.5),
		xlab = 'Minority Chamber Seat Share',
		ylab = 'Minority Committee Share',
		main = '',
		axes = FALSE)
	axis(1, at = seq(0, 0.5, 0.1))
	axis(2)
	y <- -0.05 + x
	lines(x, x, lwd = 2, lty = 2)
	lines(x, y, lwd = 6, col = rgb(0, 0, 0, 0.5))
	text(0.05, 0.45, substitute(paste(alpha, ' < 0')))
	text(0.05, 0.4, substitute(paste(beta, ' = 1')))
	legend('bottomright', legend = c('Predicted Minority Share'), lwd = 6, col = rgb(0, 0, 0, 0.5), bty = 'n')
	mtext('Majority Uniformly Overcompensated', font = 2)

	plot(1, 1, type = 'n',
		xlim = c(0, 0.5),
		ylim = c(0, 0.5),
		xlab = 'Minority Chamber Seat Share',
		ylab = 'Minority Committee Share',
		main = '',
		axes = FALSE)
	axis(1, at = seq(0, 0.5, 0.1))
	axis(2)
	y <- -0.05 + 0.7*x
	lines(x, x, lwd = 2, lty = 2)
	lines(x, y, lwd = 6, col = rgb(0, 0, 0, 0.5))
	legend('bottomright', legend = c('Predicted Minority Share'), lwd = 6, col = rgb(0, 0, 0, 0.5), bty = 'n')
	text(0.05, 0.45, substitute(paste(alpha, ' < 0')))
	text(0.05, 0.4, substitute(paste('0 < ', beta, ' < 1')))
	mtext('Majority Increasingly Overcompensated', font = 2)

	plot(1, 1, type = 'n',
		xlim = c(0, 0.5),
		ylim = c(0, 0.5),
		xlab = 'Minority Chamber Seat Share',
		ylab = 'Minority Committee Share',
		main = '',
		axes = FALSE)
	axis(1, at = seq(0, 0.5, 0.1))
	axis(2)
	y <- 0.1 + 0.7*x
	lines(x, x, lwd = 2, lty = 2)
	lines(x, y, lwd = 6, col = rgb(0, 0, 0, 0.5))
	text(0.05, 0.45, substitute(paste(alpha, ' > 0')))
	text(0.05, 0.4, substitute(paste('0 < ', beta, ' < 1')))
	legend('bottomright', legend = c('Predicted Minority Share'), lwd = 6, col = rgb(0, 0, 0, 0.5), bty = 'n')
	mtext('Majority Exploits Minority', font = 2)
	mtext('Models of Committee Seat Allocation', outer = TRUE, line = -2, cex = 1.5, font = 2)
dev.off()


## descriptive data picture (Figure 2)
	pdf('raw.pdf')
	plot(1, 1, type = 'n',
		xlim = c(0, 0.5),
		ylim = c(0, 0.5),
		main = 'Observed Relationship Between Chamber Seats and Committee Share',
		xlab = 'Minority Chamber Seat Share',
		ylab = 'Minority Committee Seat Share',
		axes = FALSE)
	axis(1)
	axis(2)
	points(jitter(data$minorityChamberShare, amount = 0.01), jitter(data$minorityCommitteeShare, amount = 0.01), col = rgb(0, 0, 0, 0.15))
	lines(c(0, 1), c(0, 1), lwd = 2, lty = 2)
	abline(lm(data$minorityCommitteeShare ~ data$minorityChamberShare), lwd = 6, col = rgb(0, 0, 0, 0.5))
	legend('bottomright', legend = c('Observed Minority Share', 'Proportionality'), lwd = c(6, 2),
		col = c(rgb(0, 0, 0, 0.5), 'black'), bty = 'n', lty = c(1, 2))
	dev.off()

## here are the descriptive statistics (Table 1) ##
	summary(data$minorityCommitteeShare)
	sd(data$minorityCommitteeShare)
	summary(data$minorityChamberShare)
	sd(data$minorityChamberShare)
	summary(data$totalCommitteeSeats)
	sd(data$totalCommitteeSeats)
	summary(data$majorityCohesiveness)
	sd(data$majorityCohesiveness)
	summary(data$proceduralPower)
	sd(data$proceduralPower)
	summary(data$competition)
	sd(data$competition)
	summary(data$proportionalityRule)
	sd(data$proportionalityRule)
	summary(data$minorityLeaderConference)
	sd(data$minorityLeaderConference)

## this estimates the main models (Table 2) ##
 	modelFe <- lm(minorityCommitteeShare ~ (minorityChamberShare)*(proportionalityRule + minorityLeaderConference + competition) + majorityCohesiveness + proceduralPower + totalCommitteeSeats + unitId, data = data)
  summary(modelFe)

 	modelRe <- lmer(minorityCommitteeShare ~ (minorityChamberShare)*(proportionalityRule + minorityLeaderConference + competition) + majorityCohesiveness + proceduralPower + totalCommitteeSeats + (1 | unitId), data = data)
	summary(modelRe)

## pull out the parameters for hypothesis testing ##
## for the fixed effects model, we aggregate ##
## all unit fixed effects for a "grand intercept" ##

	fixef <- rmnorm(1000, coef(modelFe), as.matrix(vcov(modelFe)))
	parsFix <- fixef[, 1] + rowMeans(fixef[, c(9:73)])
	parsFix <- cbind(parsFix, fixef[, c(2:8, 74:76)])
	pars <- rmnorm(1000, fixef(modelRe), as.matrix(vcov(modelRe)))
	
	withinCoef <- summary(modelFe)$coefficients[, 1][c(1:8, 74:76)]
	withinSe <- summary(modelFe)$coefficients[, 2][c(1:8, 74:76)]

	acrossCoef <- summary(modelRe)$coefficients[, 1]
	acrossSe <- summary(modelRe)$coefficients[, 2]

## this is the alpha estimate given in Table 2 for the the within model ##
	withinCoef[1] <- mean(parsFix[, 1])
	withinSe[1] <- sd(parsFix[, 1])

## this bit aggregates the critical statistics in Table 2 labeled 'Test' ##
## first the fixed effects model ##
	summary(modelFe)
	b1 <- length(parsFix[, 1][parsFix[, 1] > 0 & parsFix[, 1] < 0.5]) / length(parsFix[, 1])
	b2 <- length(parsFix[, 2][parsFix[, 2] > 0 & parsFix[, 2] < 0.5]) / length(parsFix[, 2])
	b3 <- length(parsFix[, 3][parsFix[, 3] < 0 & parsFix[, 3] > -parsFix[, 1]]) / length(parsFix[, 3])
	b4 <- length(parsFix[, 4][parsFix[, 4] < 0 & parsFix[, 4] > -parsFix[, 1]]) / length(parsFix[, 4])
	b5 <- length(parsFix[, 5][(parsFix[, 5]*max(data$competition)) < 0 & (parsFix[, 5]*max(data$competition)) > -parsFix[, 1]]) / length(parsFix[, 5])
	b6 <- length(parsFix[, 6][parsFix[, 6] > 0 & parsFix[, 6] < 0.5]) / length(parsFix[, 6])
	b7 <- length(parsFix[, 7][parsFix[, 7] > 0 & parsFix[, 7]*max(data$proceduralPower) < 0.5]) / length(parsFix[, 7])
	b8 <- length(parsFix[, 8][parsFix[, 8] > 0 & parsFix[, 8] < 0.5]) / length(parsFix[, 8])
	b9 <- length(parsFix[, 9][parsFix[, 9] > 0 & parsFix[, 9]+parsFix[, 2] < 1]) / length(parsFix[, 9])
	b10 <- length(parsFix[, 10][parsFix[, 10] > 0 & parsFix[, 10]+parsFix[, 2] < 1]) / length(parsFix[, 10])
	b11 <- length(parsFix[, 11][parsFix[, 11] > 0 & parsFix[, 11]*max(data$proceduralPower)+parsFix[, 2] < 1]) / length(parsFix[, 11])
	
	pFix <- 1 - c(b1, b2, b3, b4, b5, b6, b7, b8, b9, b10, b11)

## and now the random effects model ##
	summary(modelRe)
	b1 <- length(pars[, 1][pars[, 1] > 0 & pars[, 1] < 0.5]) / length(pars[, 1])
	b2 <- length(pars[, 2][pars[, 2] > 0 & pars[, 2] < 0.5]) / length(pars[, 2])
	b3 <- length(pars[, 3][pars[, 3] < 0 & pars[, 3] > -pars[, 1]]) / length(pars[, 3])
	b4 <- length(pars[, 4][pars[, 4] < 0 & pars[, 4] > -pars[, 1]]) / length(pars[, 4])
	b5 <- length(pars[, 5][(pars[, 5]*max(data$competition)) < 0 & (pars[, 5]*max(data$competition)) > -pars[, 1]]) / length(pars[, 5])
	b6 <- length(pars[, 6][pars[, 6] > 0 & pars[, 6] < 0.5]) / length(pars[, 6])
	b7 <- length(pars[, 7][pars[, 7] > 0 & pars[, 7]*max(data$proceduralPower) < 0.5]) / length(pars[, 7])
	b8 <- length(pars[, 8][pars[, 8] > 0 & pars[, 8] < 0.5]) / length(pars[, 8])
	b9 <- length(pars[, 9][pars[, 9] > 0 & pars[, 9]+pars[, 2] < 1]) / length(pars[, 9])
	b10 <- length(pars[, 10][pars[, 10] > 0 & pars[, 10]+pars[, 2] < 1]) / length(pars[, 10])
	b11 <- length(pars[, 11][pars[, 11] > 0 & pars[, 11]*max(data$proceduralPower)+pars[, 2] < 1]) / length(pars[, 11])
	
	p <- 1 - c(b1, b2, b3, b4, b5, b6, b7, b8, b9, b10, b11)

## these are pars, ses, and the 'test' statistics given in Table 2 ##
## note that the order is not a perfect match ##
	cbind(round(withinCoef, 3), round(withinSe, 3), round(pFix, 3), 
		round(acrossCoef, 3), round(acrossSe, 3), round(p, 3))

## now we'll create the estimated effect sizes for Figure 3 ##
## pull out the unique minority seat shares for x axis ##

	chamberSeats <- sort(unique(data$minorityChamberShare))

## we will start with competition ##
	med <- c()
	hi <- c()
	lo <- c()
	
	for(i in 1:length(chamberSeats)){
		seats <- pars[, 1] + 
			pars[, 2]* chamberSeats[i] + 
			pars[, 3]* mean(data$proportionalityRule, na.rm = TRUE) + 
			pars[, 4]* mean(data$minorityLeaderConference, na.rm = TRUE) + 
			pars[, 5]* 0.05 + 
			pars[, 6]* mean(data$majorityCohesiveness, na.rm = TRUE) + 
			pars[, 7]* mean(data$proceduralPower, na.rm = TRUE) + 
			pars[, 8]* log(mean(data$totalCommitteeSeats, na.rm = TRUE)) + 
			pars[, 9]* mean(data$proportionalityRule, na.rm = TRUE) * chamberSeats[i] + 
			pars[, 10]* mean(data$minorityLeaderConference, na.rm = TRUE)  * chamberSeats[i] + 
			pars[, 11]* 0.05 * chamberSeats[i]
			

		seatsP <- pars[, 1] + 
			pars[, 2]* chamberSeats[i] + 
			pars[, 3]* mean(data$proportionalityRule, na.rm = TRUE) + 
			pars[, 4]* mean(data$minorityLeaderConference, na.rm = TRUE) + 
			pars[, 5]* 0.25 + 
			pars[, 6]* mean(data$majorityCohesiveness, na.rm = TRUE) + 
			pars[, 7]* mean(data$proceduralPower, na.rm = TRUE) + 
			pars[, 8]* log(mean(data$totalCommitteeSeats, na.rm = TRUE)) + 
			pars[, 9]* mean(data$proportionalityRule, na.rm = TRUE) * chamberSeats[i] + 
			pars[, 10]* mean(data$minorityLeaderConference, na.rm = TRUE)  * chamberSeats[i] + 
			pars[, 11]* 0.25 * chamberSeats[i]
			
		med[i] <- quantile(seatsP - seats, 0.5)
		hi[i] <- quantile(seatsP - seats, 0.95)
		lo[i] <- quantile(seatsP - seats, 0.05)

	}			

	pdf('substantiveEffects.pdf', height = 10, width = 10)
	par(mfrow = c(2, 2))
	plot(1, 1, type = 'n',
		xlim = c(0, 0.5),
		ylim = c(-0.1, 0.1),
		ylab = 'Change in Minority Committee Share',
		xlab = 'Minority Chamber Share',
		main = 'Competition (0.05 -> 0.25)',
		axes = FALSE)
	axis(1)
	axis(2)
	lines(c(0, 0.5), c(0, 0), lty = 2)
	points(chamberSeats, med, col = rgb(0, 0, 0, 0.5), pch = 1)
	polygon(c(chamberSeats, rev(chamberSeats)), c(hi, rev(lo)), col = rgb(0, 0, 0, 0.2), border = FALSE)

## next the proportionality rule ##
	med <- c()
	hi <- c()
	lo <- c()
	
	for(i in 1:length(chamberSeats)){
		seats <- pars[, 1] + 
			pars[, 2]* chamberSeats[i] + 
			pars[, 3]* 0 + 
			pars[, 4]* mean(data$minorityLeaderConference, na.rm = TRUE) + 
			pars[, 5]* mean(data$competition, na.rm = TRUE) + 
			pars[, 6]* mean(data$majorityCohesiveness, na.rm = TRUE) + 
			pars[, 7]* mean(data$proceduralPower, na.rm = TRUE) + 
			pars[, 8]* log(mean(data$totalCommitteeSeats, na.rm = TRUE)) + 
			pars[, 9]* 0 + 
			pars[, 10]* mean(data$minorityLeaderConference, na.rm = TRUE) * chamberSeats[i] + 
			pars[, 11]* mean(data$competition, na.rm = TRUE) * chamberSeats[i]
			

		seatsP <- pars[, 1] + 
			pars[, 2]* chamberSeats[i] + 
			pars[, 3]* 1 + 
			pars[, 4]* mean(data$minorityLeaderConference, na.rm = TRUE) + 
			pars[, 5]* mean(data$competition, na.rm = TRUE) + 
			pars[, 6]* mean(data$majorityCohesiveness, na.rm = TRUE) + 
			pars[, 7]* mean(data$proceduralPower, na.rm = TRUE) + 
			pars[, 8]* log(mean(data$totalCommitteeSeats, na.rm = TRUE)) + 
			pars[, 9]* chamberSeats[i] + 
			pars[, 10]* mean(data$minorityLeaderConference, na.rm = TRUE) * chamberSeats[i] + 
			pars[, 11]* mean(data$competition, na.rm = TRUE) * chamberSeats[i]
			
		med[i] <- quantile(seatsP - seats, 0.5)
		hi[i] <- quantile(seatsP - seats, 0.95)
		lo[i] <- quantile(seatsP - seats, 0.05)

	}			

	plot(1, 1, type = 'n',
		xlim = c(0, 0.5),
		ylim = c(-0.1, 0.1),
		ylab = 'Change in Minority Committee Share',
		xlab = 'Minority Chamber Share',
		main = 'Proportionality Rule (OFF -> ON)',
		axes = FALSE)
	axis(1)
	axis(2)
	lines(c(0, 0.5), c(0, 0), lty = 2)
	points(chamberSeats, med, col = rgb(0, 0, 0, 0.5), pch = 1)
	polygon(c(chamberSeats, rev(chamberSeats)), c(hi, rev(lo)), col = rgb(0, 0, 0, 0.2), border = FALSE)

## now the conference rule ##
	med <- c()
	hi <- c()
	lo <- c()
	
	for(i in 1:length(chamberSeats)){
		seats <- pars[, 1] + 
			pars[, 2]* chamberSeats[i] + 
			pars[, 3]* mean(data$proportionalityRule, na.rm = TRUE) + 
			pars[, 4]* 0 + 
			pars[, 5]* mean(data$competition, na.rm = TRUE) + 
			pars[, 6]* mean(data$majorityCohesiveness, na.rm = TRUE) + 
			pars[, 7]* mean(data$proceduralPower, na.rm = TRUE) + 
			pars[, 8]* log(mean(data$totalCommitteeSeats, na.rm = TRUE)) + 
			pars[, 9]* mean(data$proportionalityRule, na.rm = TRUE) * chamberSeats[i] + 
			pars[, 10]* 0 + 
			pars[, 11]* mean(data$competition, na.rm = TRUE) * chamberSeats[i]
			

		seatsP <- pars[, 1] + 
			pars[, 2]* chamberSeats[i] + 
			pars[, 3]* mean(data$proportionalityRule, na.rm = TRUE) + 
			pars[, 4]* mean(data$minorityLeaderConference, na.rm = TRUE) + 
			pars[, 5]* mean(data$competition, na.rm = TRUE) + 
			pars[, 6]* mean(data$majorityCohesiveness, na.rm = TRUE) + 
			pars[, 7]* mean(data$proceduralPower, na.rm = TRUE) + 
			pars[, 8]* log(mean(data$totalCommitteeSeats, na.rm = TRUE)) + 
			pars[, 9]* mean(data$proportionalityRule, na.rm = TRUE) * chamberSeats[i] + 
			pars[, 10]* chamberSeats[i] + 
			pars[, 11]* mean(data$competition, na.rm = TRUE) * chamberSeats[i]
			
		med[i] <- quantile(seatsP - seats, 0.5)
		hi[i] <- quantile(seatsP - seats, 0.95)
		lo[i] <- quantile(seatsP - seats, 0.05)

	}			

	plot(1, 1, type = 'n',
		xlim = c(0, 0.5),
		ylim = c(-0.1, 0.1),
		ylab = 'Change in Minority Committee Share',
		xlab = 'Minority Chamber Share',
		main = 'Conference Procedure (OFF -> ON)',
		axes = FALSE)
	axis(1)
	axis(2)
	lines(c(0, 0.5), c(0, 0), lty = 2)
	points(chamberSeats, med, col = rgb(0, 0, 0, 0.5), pch = 1)
	polygon(c(chamberSeats, rev(chamberSeats)), c(hi, rev(lo)), col = rgb(0, 0, 0, 0.2), border = FALSE)


## and finally all three rules ##

	med <- c()
	hi <- c()
	lo <- c()
	medP <- c()
	hiP <- c()
	loP <- c()
	
	for(i in 1:length(chamberSeats)){
		seats <- pars[, 1] + 
			pars[, 2]* chamberSeats[i] + 
			pars[, 3]* 0 + 
			pars[, 4]* 0 + 
			pars[, 5]* 0.05 + 
			pars[, 6]* mean(data$majorityCohesiveness, na.rm = TRUE) + 
			pars[, 7]* mean(data$proceduralPower, na.rm = TRUE) + 
			pars[, 8]* log(mean(data$totalCommitteeSeats, na.rm = TRUE)) + 
			pars[, 9]* 0 + 
			pars[, 10]* 0 + 
			pars[, 11]* 0.05 * chamberSeats[i]
			

		seatsP <- pars[, 1] + 
			pars[, 2]* chamberSeats[i] + 
			pars[, 3]* 1 + 
			pars[, 4]* 1 + 
			pars[, 5]* 0.25 + 
			pars[, 6]* mean(data$majorityCohesiveness, na.rm = TRUE) + 
			pars[, 7]* mean(data$proceduralPower, na.rm = TRUE) + 
			pars[, 8]* log(mean(data$totalCommitteeSeats, na.rm = TRUE)) + 
			pars[, 9] * chamberSeats[i] + 
			pars[, 10]* chamberSeats[i] + 
			pars[, 11]* 0.25 * chamberSeats[i]
			
		med[i] <- quantile(seats, 0.5)
		hi[i] <- quantile(seats, 0.95)
		lo[i] <- quantile(seats, 0.05)
			
		medP[i] <- quantile(seatsP, 0.5)
		hiP[i] <- quantile(seatsP, 0.95)
		loP[i] <- quantile(seatsP, 0.05)

	}			


	plot(1, 1, type = 'n',
		xlim = c(0, 0.5),
		ylim = c(0, 0.5),
		ylab = 'Minority Committee Share',
		xlab = 'Minority Chamber Share',
		main = 'Best and Worst Proportionality Scenarios',
		axes = FALSE)
	axis(1)
	axis(2)
	lines(c(0, 0.5), c(0, 0.5), lty = 2)
	points(chamberSeats, med, col = rgb(0, 0, 0, 0.2), pch = 1)
	polygon(c(chamberSeats, rev(chamberSeats)), c(hi, rev(lo)), col = rgb(0, 0, 0, 0.2), border = FALSE)
	points(chamberSeats, medP, col = rgb(0, 0, 0, 0.5), pch = 3)
	polygon(c(chamberSeats, rev(chamberSeats)), c(hiP, rev(loP)), col = rgb(0, 0, 0, 0.5), border = FALSE)
	legend('bottomright', legend = c('All moderators ON', 'All moderators OFF'), border = FALSE,
		col = c(rgb(0, 0, 0, 0.5), rgb(0, 0, 0, 0.2)), bty = 'n', fill = c(rgb(0, 0, 0, 0.5), rgb(0, 0, 0, 0.2)))
	dev.off()